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Abstract. In a bipartite system subject to decoherence from two separate reservoirs, 
the entanglement is typically destroyed faster than for single reservoirs. Surprisingly 
however, the existence of separate reservoirs can also have a beneficial entangling 
effect: if the qubits are coupled and driven externally by a classical field, the system 
ends up in a stationary state characterized by a finite degree of entanglement. This 
phenomenon occurs only in a certain region of the parameter space and the structure of 
the stationary state has a universal form which does not depend on the initial state or 
on the specific physical realization of the qubits. We show that the entanglement thus 
generated can be propagated within a quantum network using simple local unitary 
operations. We suggest the use of such systems as "batteries of entanglement" in 
quantum circuits. 
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1. Introduction 

When two qubits are exposed to uncorrelated weak noise, qualitatively new phenomena, 
known as sudden death of entanglement, have been predicted theoretically [H [2] and 
confirmed experimentally for photons [3] and atoms [4]. These phenomena include the 
abrupt (rather than exponential) decay of concurrence for certain initial entangled states 
[1] and the non-additivity of decoherence rates [2]. They indicate that two different 
environments have a detrimental effect on the entanglement, which gets stronger with 
the addition of external influences, as shown by the case of driven qubits [5] and finite- 
temperature baths [6j. This contrasts to the case of a single environment, where the 
constructive role of decoherence has been known for some time [7j . Other cases in which a 
single environment can in fact support the occurence of entanglement have been pointed 
out in various contexts: for example, if the qubits are allowed to exchange excitations 
via a third continuously- monitored quantum object [8J, if the qubits are driven [9], 
dipole-coupled [10], or if vacuum fluctuations (rotating terms) are not negligible [TT] . 

In this paper we show that, for coupled and driven qubits interacting with two 
different reservoirs, entanglement can be generated and maintained at long time-scales. 
The existence of this effect does not depend on the specific physical realization of 
the qubits but requires the presence of all three ingredients: dissipation, driving, and 
coupling. We show that steady-state generation of entanglement is possible in a region 
which is outside the range of validity of the secular approximation [12], where the 
concurrence can reach a maximal universal value (half the inverse of the golden ratio) 
if a certain simple relation between the driving field and the coupling is established. 
In a wider quantum- information context, this effect could find applications such as 
on-demand creation of complex entangled states by manipulating the dissipation [13], a 
process which could be regarded as a form of adiabatic quantum computing [H] . Here we 
suggest that, since stationary entanglement is robust under decoherence, independent on 
the initial state of the qubits, relatively easy to generate, and has a system-independent 
value, such systems could be used as "entanglement batteries", and we show that, 
using simple local unitary operations available now for systems such as superconducting 
qubits, this entanglement can be harvested and transmitted further in a quantum circuit. 

2. Generation of entanglement 

We consider two qubits of Larmor frequency uj (j = 1,2), coupled by a dipole-dipole 
interaction of strength u^^{t) ^ Uj, and irradiated by external monochromatic fields 
of frequency ojj and vacuum Rabi frequencies flj. The two states of each qubit are 
denoted by |0) and |1), and for the Bell basis of the two qubits we will use the notation 
= 1/^2(101) ± |10)), and |<l>±) = l/\/2(|00) ± |11)). The dipole-dipole coupling 
between the qubits can be either direct or resulting from virtual excitations through 
a third object, as it is the case in many qubit architectures [15]. In the Schrodinger 
picture and with ^ = 1, 
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His) = ^ + co^^{t)a^,a^, + cos{co,t)a^, (1) 

i=l,2 ^ j=l,2 

where, to account for the situation in which the qubits have very different Larmor 
frequencies [16], we have considered a modulated couphng u^^it) = 2uj^^ cos[(ci;2 — uJi)t\. 
The dissipation is described by the standard Born-Markov master equation [12j . 

pis) = -i[H^^)^p^^)]+C[p^% (2) 

where the Liouvillean 

^[p^^^] = E ^ (2a7p(^v; - a>7p(^) - p(^)a;.-) (3) 

i=i,2 ^ 

models the longitudinal dampings of the qubits, and T j is the standard energy relaxation 
rate. 

We work in a rotating reference frame, characterized by the transformation R = 
exp[i{ujiaf + U20'2)t/2]. Using the identities K^afR = af exp{^iujjt) and eliminating 
the fast counter-rotating terms via the rotating wave approximation, we obtain a time- 
independent master equation in the rotating frame, 

p = -t[H,p]+C[p], (4) 

with p = Rp^^^R\ and 

H - E (1^1 + ) + ^ (^?^2 + ^r^!) • (5) 

Here Sj = Vj — Uj are the detunings of the qubits from the corresponding driving 
frequencies. In order to use the rotating wave approximation, we have taken ui,U2, \uji — 
UJ2\ ^ Qj,6j,uj^^,rj. For the case of qubits with close Larmor frequencies, there is no 
need to modulate the coupling: an effective Hamiltonian of the type Eq. (|5]) can be 
obtained [T7] by the same transformations and working at equal driving frequencies 
Ui = bJ2- From now on, we will also refer to the interaction part (containing w^^) of the 
Hamiltonian Eq. ^ as if^^, and the rest as [H = Hq + H^^). 

To explore the entanglement properties of this system, we have first solved 
numerically Eq. (jl]) for 6j = 6, Qj = Q, and Tj = F. The entanglement between the 
qubits at any time is characterized by Woofers' concurrence [H], defined as max{0,C}, 
where C = Ai — A2 — A3 — A4, and AjS are the eigenvalues of ^Jy/ppy/p in decreasing order, 
with p = {(jy ® o"^)p*(o"^ ® a^). We find that for most of the values of fl and u^^ either 
sudden death of entanglement or exponential decay of entanglement occurs, depending 
on the initial state [1]. However, for certain values of the driving field fl and coupling 
^xx ^ process of concurrence buildup occurs, as shown in Fig. [1] for the case of initial 
Werner states [H], defined as pw = [(1 - /)/3]/+ [(4/ - l)/3]|^-)(^-|. 

We first note that this effect is qualitatively different from the case of a single 
reservoir, with Liouvillean superoperator 

Cs^ngle[p] = ^ {2S- pS+ - S+ p - pS+ S') , (6) 
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with S*^ = J2j=i,2 being collective spin operators. In this situation, the emergence 
of an entangled steady state via evolution depends on the initial state, as shown in Fig. 
121 where we plotted the time-dependent concurrence for three classes of states, Werner 
states, Yu-Eberly (YE) states [22] (states of the type pye = (2/3)|^+)(^+| + [(1 - 
a)/3]|ll)(ll| + (a/3)|00))00|) , and states diagonal in the subspace with one excitation, 
which we denote by Peg-ge = (1 ~ a)|10)(10| + a|01)(01|. The reason for this is the 
existence, for common reservoirs, of a decoherence-free subspace: as a result, the amount 
of entanglement contained in the asymptotic state depends in general on the projection 
of the initial state onto this subspace [20]. A number of effects related to decoherence- 
free subspaces in interacting qubits have been studied recently: for example, it has been 
shown how to use the qubit-qubit coupling to effectively produce a low-decoherence 
subspace even when the qubits have additional separate decoherence channels [21], and 
how to create maximally-entangled state by suitable measurements on the decay photons 




Figure 2. (color online). Concurrence evolution for (left) YE states, for (middle) 
Werner states, and (right) for eg — ge states. For all the plots, w^^ = lOF, = l.SF. 



We now return to the case of two reservoirs. Entanglement generation occurs only 
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in a certain region of the parameter space: in Fig. [3] we have plotted the concurrence 
for a long time-scale evolution t ^ r2~^,r~^. The maximum value of the concurrence 
is reached at resonance {6 = 0) and it increases slowly with Q/T and u^^/T. Although 
intuitively one might expect that increasing either the coupling or the pumping strength 
would increase the steady-state concurrence, this is not the case: both the ratios Q/T 
and u^^/T are important, showing that all three processes (coupling, pumping, and 
decoherence) contribute to this effect. It also means that the secular approximation [12] 
is not valid in the region of the parameter space where the effect occurs. 

To get some insight into the mechanism responsible for the generation of steady- 
state concurrence, we will take a perturbative approach: if there is no coupling between 
the two qubits, then in the stationary regime the bipartite system is described by 
the separable density matrix Pi°°''p2°°''- Consider a single qubit j under constant on- 
resonance driving: in the rotating frame, the master equation is simply 

Pj = - yt^^ Pj] + 2 (2^"Pi^^ - fx+ff-pj - PifT+o--) , (7) 



and the steady-state is (see [Appendix A] ) 

joo) _ 1 ( \ 

Pj r2 + 1^ -inr r^ + n^ )' ^ ^ 

We now aim at finding the first-order contribution in u^^/T ^ 1 to the density 
matrix and the concurrence. We write the solution as p*^°°-' ~ Pi^^P2^^ + P^^, where, in 
order to ensure a unit value for p^°°^ , we search for a traceless density matrix p^^ which 
satisfies the equation —i[Ho,p^^] — i[H^^ , p^^^ p''^^] + £p^^ = 0. Here we neglect the 
second-order contribution in u^^/T coming from the commutator [H^^,p^^]. With the 
notation p^^,im = {jk\p^^\l''TT'), we find that p^^ has matrix elements Pii,oo = Poo,ii ^ 
2zfi2ra;-/(r2 + 2ny, = p-*o = p-oo = p-*i = 2nT'u^yiT' + 2V?f and 

the rest zero. The quantity C associated to the density matrix Pi°°V2°°^ + P^^ is 
C = 2r]2(2a;^^r - n^){T^ + 2n^)-^, and we have checked that it gives excellent fits 
to the numerical results in the region where u^^/T ^ 1 for both negative and positive 
values of C; it also predicts that C will become positive for u^^ = Q'^/2r, in agreement 
with the exact result derived below Eq. ffTOj) . 

Although in the region of interest (finite, large concurrence) the problem is clearly 
nonperturbative, the perturbative solution serves as a heuristic guide to find the 
stationary density matrix that solves Eq. (jlj) under the condition p = 0. Note that this 
condition implies solving for the 15 unkown real entries of p = 0. Surprisingly, although 
the algebraic calculations are rather complicated, at resonance we find a simple solution 
for p(°°), 

Joo) ^ (r2 + 2ny r n('^)n(^) + n^A 

^ [(r2 + 2fi2)2+4^a;a;2r2]2 V^l P^ ^P ) 

p2, ,xx2 

+ — (1 - (T«)(l - ^(2)). 

[(r2 + 21]2)2+4^xx2r2]2^ A ^) 

This solution preserves the structure of the perturbative case - note that the role of the 
second term is just to ensure that Trp^°°^ = 1, due to the appearance of the prefactor 
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Figure 3. (color online). Numerical values for the concurrence at a time t — lO'^F^^, 
as a function of 57 and o;^^. The dashed line is a plot of the parabola Eq. (|10p and the 
dotted line is a plot of Eq (|12p . 



in front of Pi°°''p2°°'' + P^^ the first term. Similarly, C can be calculated analytically 
and the result resembles the nonperturbative case, 

(r2 + 2fi2)2 + 4u;-2r2- 

The crossover between the region of finite steady-state concurrence and the region 
of zero concurrence is given by the parabola 

and the concurrence reaches a maximum value of 



(10) 



c 



for 



Q'^ 1 / 

2r 2r ^ ^ ^ 
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(12) 



At large values of the driving field Q/T ^ 1, we obtain a maximum concurrence of 
1/(1 + a/5) = 0.309 (which incidentally is half of the inverse of the golden ratio) , a 
result which we have confirmed numerically; the density matrix corresponding to this 
concurrence is 

/ i \ 



p(oo) 
r max 



1 



2y5(2v^- 



-/ + 



1 



2v^ 





-i 1 + J 



(13) 
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In the case of finite detuning, the calculations are more involved but it is still 
possible to obtain analytical expressions. We find for C 

C = {2u-^\r\ - Q^) , (14) 

where 

AT = [(|f|2 + 2^]2)2 + 4a;--|f|2(a;-- + 25)|]^/^ (15) 

and r = r + 2i6. The fidelities with respect to the Bell basis 
also have simple forms F(|V'+), p(°°)) = QAf-^ ^2\t\^ + , F(|V'"), p^^^) = 
^W-^ F(|0+),p(°°)) = ^1/2 - l]W-4(2r2 + A6u^^), and F(|0-), p(°°)) = 

y^l/2 - fiW-4(2r2 + fi2 + 45a;^^). In Fig. H we plot the fidelities at resonance for 
various fi's. 



F 




2 4 6 8 10 

Figure 4. The fidelities F{\il;+) , p'^°°'>) (gray), (red), F(|0+), p(°°)) 

(blue), and i^(|0^), p(°°)) (green), for = 7r and 6^0. 

The analysis above shows that the mechanism of generating entanglement is related 
to the existence of off-diagonal matrix elements (single-qubit coherences) in p^^^'- these 
get coupled by the interaction (as shown by the existence of linear terms in F and 
fl in p^^). It is interesting to note that the role of driving is solely to pump energy 
in the system: the source can be a classical one and perfect coherence is not a strong 
requirement either. Indeed, if the qubits are driven by incoherent sources of components 
k, J2k ^'■'^^ cos{ijjjt + (p^''^){af + 0"! ), then the same formula for the concurrence Eq. (1141) 
can be obtained, with = \J2k^^''^ ^wi^V^''^)]^ and, as long as this quantity is not 
exactly zero, the effect described above will be obtained. 

3. Propagation of entanglement in quantum circuits 

We now show that designing quantum circuits that would harvest and propagate 
this stationary entanglement at places where it is needed in a quantum processor is 
possible using a tunable Jaynes-Cummings interaction, a technology already available 
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for example in the field of superconducting qubits [231 El] • Such schemes are of general 
interest for entanglement distribution between the nodes of a quantum networks |25] . 
Since entanglement can be regarded as a resource for processing information, such a 
circuit could effectively work as an "entanglement battery": when entangled fields are 
needed in other parts of the processor, it would be enough to outcouple them using for 
example two transmission lines, similar to standard voltage or current sources, as shown 
schematically in Fig. |5l 




Figure 5. Transmission of entanglement: schematic of the two coupled qubits and 
the cavities. 

Based on the already-achieved qubit manipulation protocols in these quantum 
computing architectures, we consider a circuit design in which the qubits are coupled 
to each other and also to two sections of two coplanar waveguides (or striplines) used 
as resonating cavities. Then, after the system reaches the stationary state, and for a 
short time r ^ u^^~^,T~^,Q~^ , we couple simultaneously the qubits 1,2 to the modes 
a and b, respectively, of the two cavities. Alternatively, a tunable interqubit coupling 
can be used to turn off u^^ [26]. The coupling strength (vacuum Rabi frequency) g 
of each qubit-resonator system is chosen such that g ^ u^^, F, f2. Then the effective 
Hamiltonian describing the oscillations between a qubit and its corresponding cavity is, 
in the rotating frame, of the Jaynes-Cummings type {g{a~^ai +a~af) for the first qubit, 
and g{b^o-2 +h~a2) for the second one). 

The dynamic evolution is, under the conditions specified above, such that vacuum 
Rabi oscillation between the qubits and the respective resonators occur, which can be 



described in a simple way as follows: 

|li,Oa) ^ cos(^r)|li,Oa) -isin(^r)|Oi,la), (16) 

|0i,0,) ^ |0i,0,), (17) 

|l2,0b) ^ cos((?r)|l2,0b) -^sin((7r)|02,lfe), (18) 

|02,0b) ^ |02,0fe), (19) 



where for clarity we have introduce the subscripts 1 and 2 to denote the two qubits. We 
see that for gr = tt /2 this transformation realizes two independent 7r/2 rotations in the 



Generation and propagation of entanglement 



9 



(21) 



two-dimensional subspaces spanned by |Oi,la), |li,Oa) and |02,lf,), |l2,0fe) respectively. 
As a result, a maximally entangled state for example 1^1'"*') will be transferred, after a 
time r, into a photonic Bell state 

|^+) ® \OaOb) ^ -«|0l02) ® ^(llaOfe) + |0J,)), (20) 

where |la), and |0a) are the two states of the cavity a, and similar for b. In general, for 
an arbitrary r, the density matrix in the photon basis {\la, U), |la, Of,), \0a, lb), |0a, Ob)} 
can be obtained by applying the transformation Eq. ([16]- [H]) and tracing out the qubits' 
degrees of freedom. The result is 

/ \ 

(l/2)sin2^r (l/2)cos2^r 

(l/2)sin2^r (l/2)cos2^r 

V cos^^r J 

with concurrence C given by C = sin^ gr. 

We now apply the same procedure for the general case in which the initial 
density matrix is the steady-state solution for two interacting qubits and two 
independent reservoirs. If we denote the elements of this matrix by pfkj^"'^^^^'^^ = 
{jakb\p^^'^^^^^^^^^'^^\la,'m'b) , where j,k,l,m G {0,1}, we find, after a qubit-resonator 
coupling time r and after tracing out the qubits, the following density matrix elements 
(we give only the values of the 9 independent elements, the rest can be found 
from 2^7- (p (transmitted)^ ^ ]^ Pjfc,!m™'"'''^^ = ~ Pb^'^jk^'^^'"^^* ) ^hc photou basis 
{|l„lb),|l„Ofe),|0„lfe),|0„0,)}, 

p£ir'"^'^ = (sin^7r)Vi5i (22) 

^(t„ittcd)^ -^(sin^r)VST:;o (23) 

^(t— ed)^ -^(sin^,r)VSSi (24) 

^(t— ed)^ -(sin^r)VSo (25) 

^(tmnsmitted) ^ ^^.^ C/t) Vi^O + (sm ^r)2(cOS ^t) ^1^1 (26) 

pSor'*^'^ = (sin^7r)Vai (27) 
^(tjansmitted) ^ _ ^^^.^ c/r)p^~|Jo - i(sin gr) (cos ^r) (28) 

^(^t^ransmitted) ^ ^^.^ ^t) ^0^01 + (sin grficOS ^t) VlT^U (29) 
(transmitted) ■/ ■ \2 (oo) •/ . \/ \2 (oo) /on\ 

Poi.oo = -Hsm^r) p^i ,^o-«(sm ^r) (cos ^r) jo (30) 

The concurrence corresponding to this density matrix is plotted in Fig. |6] for 
various qubit-resonator coupling times. We also note that an interesting property of 
the procedure described above is that, if the input is a density matrix in the X-form 
(nonzero elements only on the diagonal and counterdiagonal [UEIEIIS]), the transmitted 
density matrix is also in the X-form. 

We find, after a qubit-resonator coupling time gr = 7r/2, and after tracing out the 
qubits, the following density matrix in the photon basis {11^1^), |laOb), |Oalft), |OaOfe)}, 

transmitted _ / '\\j+k-j+k+l+m {'X>) /Ql N 

Pjk,lm — \ ^) '■ Pjk,lm- > 
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Figure 6. Photon-photon concurrence as a function of couphng time r for lo^^ — 20T 
and (5 = 0. 

This density matrix has the same concurrence as the initial one, in fact, it is identical to 
the initial qubit-qubit density matrix up to the local transformations \0a) — |0a), |la) — ^ 
—i\la) and similarly for photon b. These are phase quantum gates, and - if the exact 
form of the steady-state density matrix needs to be recovered - can be implemented 
by a global redefinition of the phase and two separate vr rotations around the z-axis, 
exp[-i7c] exp[-i7r(cr^ + cr^)/2]. 

Repeating this procedure with further circuit elements, entangled pairs can be 
outcoupled and processed in various parts of a quantum processor [27] , where they can 
be distilled [28] to maximally entangled states (used for example in quantum gates) or 
they can be used directly for tasks that do not require highly entangled states, such 
as quantum catalysis, - making possible certain local transformations and enhancing 
distillation algorithms [29]. Also, we suggest that since from Eq. ( IT3l) is relatively 
easy to achieve and has a universal form (independent on any qubit parameter or physical 
realization), it can also serve as a standard of entanglement for quantum networks. In 
other words, since entanglement is a measurable quantity, it should have a metrological 
standard associated: one can go further and speculate that the value of concurrence 
of 0.309 could then play the same role in metrology as other numbers do for the 
corresponding physical quantities {e.g. the temperature of the triple point of water, 
the value of the transition frequency between the two hyperfine states of Cs-133 atoms 
in atomic clocks, etc.). 

4. Conclusions 

In conclusion, we have shown that entanglement emerges in a certain region of the 
parameter space of driven, coupled qubits interacting with separate reservoirs; we give 
analytical results characterizing the stationary entanglement and we show how to further 
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propagate this entanglement in a quantum network. 
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Appendix A. Single qubit case 

Here we review the case of one qubit under dissipation. The rotating-frame single-qubit 
density matrix pj can be parametrized on the Bloch-sphere, 

Pj = lil + ra), (A.l) 

where r is the Bloch vector of components {rxi^yi^z)-, the norm of which indicates the 
purity of the state, Trp^j = (1 + |r|^)/2. The decoherence is described by Eq. ([7]), 

Pi = - yt^^ Pj\ + 2 (2^"Pi^^ - (y^(y'pj - Pi^^f^") • (A.2) 

From Eq. (IA.lfA.2p we obtain following kinetic equation for the components of the Bloch 
vector, 

r 

The last equation has a simple solution, r^it) = r^(0) exp(— rt/2). The first two 
equations can be solved by introducing the notation \C,{t)) = [rz,ry]'^ and |x) = [— T^O]; 
we find the equation 

j^\m) = M\m + \x). (A.3) 

Here the matrix M is defined by 



nry r(l -l-r^), f'y r,^y ^^2! 



and has eigenvalues Ai^2 = (1/4) (—3r =f vT^-TgTF) corresponding to (unnormalized) 
eigenvectors |^i^2) = [(r ± — 16fi^)/4i7, 1]-'". The stationary solutions of Eq. (1A.3P 
are obtained for t ^ F^^ as 1^^°°'') = —M~^\x), resulting in the density matrix Eq. 

(00) _ / r2+2Q2 r2+2f72 1 ( \ 

Pj - inr r2+!^2 . I,a.oj 
V r2+2Q2 r2+2f72 / 
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P 



(A.6) 



The elements of this density matrix are shown in Fig. IA1[ In the Bloch sphere 
representation, the stationary state density matrix is parametrized by r^°°^ = 0, 
r(°°) = 2fi^/(2^]2 + r2), and r'f^^ = -T^/iV^ + 2n^), describing a semicircle if radius 
1/2 in the y — O — z plane, (r^ + 1/2)^ + r^/2 = 1/4 (shown in blue in the Bloch sphere 

of Fig. m . 

Depending on the values of Q/T the eigenvalues Ai^2 can have imaginary parts or 
can be real, corresponding respectively to underdamped and overdamped dynamics. In 
the underdamped regime, defined by Q/T > 1/4, the elements of the density matrix 
oscillate (Rabi oscillations) before reaching the stationary state; deep in this regime 
when f2/r ^ 1/4, we find 

' 1/2 -iT/2n \ 
^ iT/2n 1/2 ) 

In the overdamped regime defined by ^7/^ < 1/4 there are no Rabi oscillations; deep in 
this regime, when Q/T ^ 1/4, we can write approximately 

' -in/T \ 

^ in/T 1 ) ■ 

Finally, the critically damped regime occurs at = r/4, and 

' 1/18 -i2/9 \ 
^ i2/9 17/18 ) ■ 

From this, we see that for values of u^^ of the order of F the emergence of stationary 
entangled states depends on having relatively large single-qubit coherences r^, which 
happens not far from the critically damped regime (the single-qubit coherence reaches 
a maximum of l/2-\/2 at fi/F = l/\/2). This is due to the fact that the interaction 
couples the qubits via the off-diagonal terms. 



P 



(<X)) 



p 



(oo) 



(A.7) 



(A. 





Figure Al. The elements of the steady state of the single qubit as a function of fl/T 
(left) and the corresponding representation (semicircle, in blue) on the Bloch sphere 
(right). 
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